# Copyright (c) 2008-2011, Thomas Handorf
# 
# This file is part of libScopes.
# libScopes is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# libScopes is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with libScopes.  If not, see <http://www.gnu.org/licenses/>.
#
package Scopes::Net;

use Scopes;
use strict;
use warnings;

=head1 Scopes::Net

Module for holding a metabolic reaction network.

=cut

use Data::Dumper;
use overload '""' => \&toString;
use constant NOOWNER => 1;
use Term::ANSIColor qw(:constants);
use Sets;
use Dmp;
sub new{
   my $class=shift;
   return bless {c_net=>new Scopes::c_net(),comID=>[],comhash=>{},reaID=>[],reahash=>{},chnhash=>{},chnID=>[],compartments=>{},comIDtoComp=>{},clabels=>{},rlabels=>{}},$class;
}
# sub STORABLE_freeze{
#    my ($self, $cloning) = @_;
#    return if $cloning;         # Regular default serialization
#    return "",@{$self}{cnet,comID,comhash,reaID,reahash,compartments,comIDtoComp};
# }
# sub STORABLE_thaw{
#    my $self=shift;
#    my $cloning=shift;
#  my $str=shift;
#  @{$self}{cnet,comID,comhash,reaID,reahash,compartments,comIDtoComp}=@_;
# }


=head2 C<addRea> (method) 

=over

=over

=item $Net->addRea($substratesHash,$productsHash,[$reversibility,[$catalysts,[inhibitors,[inhibited,[propagators]]]]])

=item $Net->addRea($substrates,$products,[$reversibility,[$stoich_subs,$stoich_prds,[$catalysts,[inhibitors,[inhibited,[propagators]]]]]])

=back

Add a reaction without specifying a reaction ID. $substrates and $products are arrays with substrate and product IDs, $reversibility is the reversibility of the reaction (indicated by -1 (backward direction), 0 (reversible=default) and 1 (forward direction), and $stoich_subs,$stoich_prds are the stoichiometric coefficients of the substrates and products given as arrays with the same order as $substrates and $products. If no stoichiometry is given, 1 is assumed for all compounds.
Note: if you supply hashes, don't supply the stoichiometry arrays

=head3 Example: 2 C00002 + C02308 <=> 2 C00008 + C02307

   $Net->addRea([C02308,C00002],[C02307,C00008],0,[1,2],[1,2]);
   $Net->addRea({C02308=>1,C00002=>2},{C02307=>1,C00008=>2},0);

=cut

sub addRea{
   my $self=shift;
   return $self->replaceRea(-1,@_);
}

=head2 C<replaceRea> (method) 

=over

=item $Net->replaceRea($reaID,$substratesHash,$productsHash,[$reversibility,[$catalysts,[inhibitors,[inhibited,[propagators]]]]])

=item $Net->replaceRea($reaID,$substrates,$products,[$reversibility,[$stoich_subs,$stoich_prds,[$catalysts,[inhibitors,[inhibited,[propagators]]]]]])

=back

Add a reaction with specifying a reaction ID. $reaID is the ID of the added reaction. $substrates and $products are arrays with substrate and product IDs, $reversibility is the reversibility of the reaction (indicated by -1 (backward direction), 0 (reversible=default) and 1 (forward direction), and $stoich_subs,$stoich_prds are the stoichiometric coefficients of the substrates and products given as arrays with the same order as $substrates and $products. If no stoichiometry is given, 1 is assumed for all compounds.
Note: if you supply hashes, don't supply the stoichiometry arrays

=head3 Example: 2 C00002 + C02308 <=> 2 C00008 + C02307

   $Net->replaceRea(R00076,[C02308,C00002],[C02307,C00008],0,[1,2],[1,2]);
   $Net->replaceRea(R00076,{C02308=>1,C00002=>2},{C02307=>1,C00008=>2},0);

=cut

sub replaceRea{
   my $self=shift;
   my @par=@_;
   my $add=0;
   my $rid=$par[0]; # given rea id
#  print Dumper(\@_);
   {
      no warnings;
      $par[0]=$self->_getreaidx($par[0]) unless $par[0]==-1; # addRea sets reaID to -1
   }
#   print "replaceRea:",Dumper(\@par);
   if (ref($par[1]) eq 'HASH'){ # substrate handling for hashes
      push(@par,0) if scalar(@par)==3; # no rev given
      splice(@par,4,0,[]); # add stoich_subs
      my $coms=[];
      for my $k (keys %{$par[1]}){
         push(@$coms,$self->_getcomidx($k));
         push(@{$par[4]},$par[1]->{$k});
      }
      $par[1]=_cIXAR($coms);
      $par[4]=_cIXARF($par[4]);
   } else {
      $par[1]=$self->_getcomar($par[1]);
   }
   if (ref($par[2]) eq 'HASH'){ # product handling for hashes
      push(@par,0) if scalar(@par)==3; # no rev given
      splice(@par,5,0,[]); # add stoich_prds
      my $coms=[];
      for my $k (keys %{$par[2]}){
         push(@$coms,$self->_getcomidx($k));
         push(@{$par[5]},$par[2]->{$k});
      }
      $par[2]=_cIXAR($coms);
      $par[5]=_cIXARF($par[5]);
   } else {
      $par[2]=$self->_getcomar($par[2]);
   }
   $par[6]=$self->_getcomar($par[6]) if $par[6];
   $par[7]=$self->_getcomar($par[7]) if $par[7];
   $par[8]=$self->_getcomar($par[8]) if $par[8];
   $par[9]=$self->_getcomar($par[9]) if $par[9];
   if ($par[0]==-1){ # addRea
      $add=1;
      shift(@par);
   }
#  print "replaceRea: $rid ",Dumper(\@par);
   my $idx=($add ? $self->{c_net}->addRea(@par): $self->{c_net}->replaceRea(@par));
   if ($idx<0){
      print STDERR "Error: reaction not added\n";
   }
   if (!$self->{reaID}->[$idx]){ # reaID not yet defined -> create an ID
      my $cc=$self->numR;
      while (exists $self->{reahash}->{"R".sprintf("%05d",$cc)}) {$cc++};
      $self->{reaID}->[$idx]="R".sprintf("%05d",$cc);
      $self->{reahash}->{"R".sprintf("%05d",$cc)}=$idx;
   }
}

=head2 C<create> (method) 

=over

=item $Net->create($reaIDs,$comIDs,$reaDefs)


=back

creates a whole network at once. $reaIDs and comIDs are arrays with reaction and compound IDs. They also determine the order (and indices) of the compounds and reactions, in particular also the order in the scope arrays. $reaDefs contains an array with definitions for each reaction. a definition is an array analogously to the paramters of replaceRea, i.e. ($reaID,$substrates,$products[,$rev[,$subs_stoich,$prods_stoich]])


=head3 Example: 2 C00002 + C02308 <=> 2 C00008 + C02307; C00002 => C00008 + C00009

   $Net->create([C00002,C00008,C00009,C02307,C02308],
                [R00076,R00022],
                   [
                     [R00076,[C02308,C00002],[C02307,C00008],0,[1,2],[1,2]],
                     [R00022,[C00002],[C00008,C00009],1,[1],[1,1]]
                   ]
                ]);

=cut


sub create{
   my $self=shift;
   my $reas=shift;
   my $coms=shift;
   my $defs=shift;
   $self->{c_net}=new Scopes::c_net();
   $self->{reaID}=$reas;
   $self->{reahash}={};
   my $rc=0;
   for my $r (@$reas){
      $self->{reahash}->{$r}=$rc++;
   }
   $self->{comID}=$coms;
   $self->{comhash}={};
   my $cc=0;
   for my $c (@$coms){
      $self->{comhash}->{$c}=$cc++;
   }
   $self->{c_net}->resize($rc,$cc);
   for my $def (@$defs){
      my ($rea,$subs,$prds,$rev,$subs_st,$prds_st)=@$def;
      $rea=$self->_getreaidx($rea);
      $subs=$self->_getcomar($subs);
      $prds=$self->_getcomar($prds);
      my $ret=-1;
      $ret=$self->{c_net}->replaceRea($rea,$subs,$prds) if scalar(@$def)==3;
      $ret=$self->{c_net}->replaceRea($rea,$subs,$prds,$rev) if scalar(@$def)==4;
      $ret=$self->{c_net}->replaceRea($rea,$subs,$prds,$rev,$subs_st,$prds_st) if scalar(@$def)==6;
      if ($ret<0){
         print STDERR "Error rea not added\n";
      }
   }
}

=head2 C<setMembers> (method) 

defines family members of family compounds.

=over

=item $Net->setMembers($comIDs)

=cut

sub setMembers{
   my $self=shift;
   my $n=shift;
   my $mems=shift;
   $mems=$self->_getcomar($mems);
   $self->{c_net}->setMembers($self->_getcomidx($n),$mems);
}

=head2 C<setChannels> (method) 

sets channels for compounds. Channels define which members of different compound families are related. An example 
would be to use uniprotids as channels. Families which represent phosphorylated and unphosphorylated versions of 
proteins can thereby relate the phosphorylate and unphosphorylated form of the members.

=over

=item $Net->setChannels($ChannelIDs)

=cut

sub setChannels{
#   print "setChannels", Dumper \@_;
   my $self=shift;
   my $n=shift;
   my $chns=shift;
   
   $chns=$self->_getchnar($chns);
   $self->{c_net}->setChannels($self->_getcomidx($n),$chns);
   
}

=head2 C<size> (method) 

returns the size (number of reactions, number of compounds) of the network

=cut

sub size{
   my $self=shift;
   return ($self->{c_net}->numR,$self->{c_net}->numC);
}

=head2 C<numR> (method) 

returns the number of reactions of the network

=cut

sub numR{
   my $self=shift;
   return $self->{c_net}->numR;
}

=head2 C<numC> (method) 

returns the number of compounds of the network

=cut

sub numC{
   my $self=shift;
   return $self->{c_net}->numC;
}

=head2 C<at> (method) 

=over

=item $Net->at($reaIdx,$comIdx$)

=back

returns the entry of the stoichiometric matrix at a position. The index of a reaction and of a compound have to be provided.

=cut

sub at{
   return $_[0]->{c_net}->at($_[1],$_[2]);
}

sub print_matrix{
   return $_[0]->{c_net}->print_matrix;
}
sub revIX{
   return $_[0]->{c_net}->rev($_[1]);
}
sub revID{
   return $_[0]->{c_net}->rev($_[0]->{reahash}->{$_[1]});
}
sub reaMask{
   return $_[0]->{c_net}->{reaMask};
}
sub selectID{
   my ($self,$id)=@_;
   return if ! exists $self->{reahash}->{$id};
   $self->{c_net}->{reaMask}->put($self->{reahash}->{$id},1);
}
sub unselectID{
   my ($self,$id)=@_;
   return if ! exists $self->{reahash}->{$id};
   $self->{c_net}->{reaMask}->put($self->{reahash}->{$id},0);
}
sub selectedID{
   my ($self,$id)=@_;
   return 0 if ! exists $self->{reahash}->{$id};
   return $self->{c_net}->{reaMask}->at($self->{reahash}->{$id});
}
sub selectedIX{
   my ($self,$ix)=@_;
   return 0 if ! $self->{reaID}->[$ix];
   return $self->{c_net}->{reaMask}->at($ix);
}

#-- networks inspection -----

sub neighborsIX{
   $_[0]->{c_net}->neighbors($_[1])->toa();
}
sub succsIX{
   return $_[0]->{c_net}->succs($_[1])->toa();
}
sub predsIX{
   $_[0]->{c_net}->preds($_[1])->toa();
}
sub metabolitesIX{
   $_[0]->{c_net}->metabolites($_[1])->toa();
}
sub productsIX{
   $_[0]->{c_net}->products($_[1])->toa();
}
sub substratesIX{
   $_[0]->{c_net}->substrates($_[1])->toa();
}
sub prod_stoichIX{
   $_[0]->{c_net}->prod_stoich($_[1])->toa();
}
sub subs_stoichIX{
   $_[0]->{c_net}->subs_stoich($_[1])->toa();
}
sub catalystsIX{
   $_[0]->{c_net}->catalysts($_[1])->toa();
}
sub inhibitorsIX{
   $_[0]->{c_net}->inhibitors($_[1])->toa();
}
sub inhibitedIX{
   $_[0]->{c_net}->inhibited($_[1])->toa();
}
sub propagatorsIX{
   $_[0]->{c_net}->propagators($_[1])->toa();
}
sub membersIX{
   $_[0]->{c_net}->members($_[1])->toa();
}
sub channelsIX{
   $_[0]->{c_net}->channels($_[1])->toa();
}
sub neighborsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{comhash}->{$id};
   return $self->rIXID($self->neighborsIX($self->{comhash}->{$id}));
}
sub succsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{comhash}->{$id};
   return $self->rIXID($self->succsIX($self->{comhash}->{$id}));
}
sub predsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{comhash}->{$id};
   return $self->rIXID($self->predsIX($self->{comhash}->{$id}));
}
sub metabolitesID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->metabolitesIX($self->{reahash}->{$id}));
}
sub productsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->productsIX($self->{reahash}->{$id}));
}
sub substratesID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->substratesIX($self->{reahash}->{$id}));
}
sub prod_stoichID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->prod_stoichIX($self->{reahash}->{$id});
}
sub subs_stoichID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->subs_stoichIX($self->{reahash}->{$id});
}
sub catalystsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->catalystsIX($self->{reahash}->{$id}));
}
sub inhibitorsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->inhibitorsIX($self->{reahash}->{$id}));
}
sub propagatorsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->propagatorsIX($self->{reahash}->{$id}));
}
sub inhibitedID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{reahash}->{$id};
   return $self->cIXID($self->inhibitedIX($self->{reahash}->{$id}));
}
sub membersID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{comhash}->{$id};
   return $self->cIXID($self->membersIX($self->{comhash}->{$id}));
}
sub channelsID{
   my $self=shift;
   my $id=shift;
   return [] if ! exists $self->{comhash}->{$id};
   return $self->chIXID($self->channelsIX($self->{comhash}->{$id}));
}

sub reaDegree{
   my $self=shift;
   my $id=shift;
   return $self->{c_net}->substrates($self->{reahash}->{$id})->size()+$self->{c_net}->products($self->{reahash}->{$id})->size();
}
sub reaInDegree{
   my $self=shift;
   my $id=shift;
   return $self->{c_net}->substrates($self->{reahash}->{$id})->size();
}
sub reaOutDegree{
   my $self=shift;
   my $id=shift;
   return $self->{c_net}->products($self->{reahash}->{$id})->size();
}
sub comDegree{
   my $self=shift;
   my $id=shift;
   return $self->{c_net}->preds($self->{comhash}->{$id})->size()+$self->{c_net}->succs($self->{comhash}->{$id})->size();
}
sub comInDegree{
   my $self=shift;
   my $id=shift;
   return $self->{c_net}->preds($self->{comhash}->{$id})->size();
}
sub comOutDegree{
   my $self=shift;
   my $id=shift;
   return $self->{c_net}->succs($self->{comhash}->{$id})->size();
}
sub reaexists{
   return 0; # not yet implemented
}

=head2 searchReas (method)

=over

=item $Net->searchReas($query)

=back

searches for all reactions fullfilling the condition in query. A condition is defined as follows:

condition: cond_list [ sidesep cond_list ]

cond_list: and_cond [ "|" cond_list]

and_cond : com_cond [ " " and_cond]

com_cond : comID | !comID

sidesep  : "<=>"|"<="|"=>"

A condition contains compounds IDs (comID) which are linked by various AND" ", OR"|" and NOT"!" operations. The separator sidesep defines whether the left condition only on the substrate side and the right condition only at the products side of a reaction (=>) , vice versa (<=) or both (<=>). If ommitted the one condition is simply checked on both sides.

Retuns and arrayref containing the reactions which met the condition.

=head3 Example:

   $Net->searchReas("C00002<=>C00008"); #all reactions with C00002 on one side and C00008 on the otherside
   $Net->searchReas("C00002=>C00008"); #dito, but C00002 must be substrate
   $Net->searchReas("C00002<=>C00008|C00046"); #all reactions with C00002 on one side and   on the otherside
   $Net->searchReas("C00002 C00001<=>C00008|C00046"); #all reactions with (C00002 and C00001) on one side and (C00008 or C00046) on the otherside

=cut

sub searchReas{
   #$net->searchRea("C00002|C00003 <=> !C00001 C00004")
   my $self=shift;
   my $q=shift; # query string
   $q=~s/^\s*//;
   $q=~s/\s*$//;
   print "search $q\n";
   my ($lside,$dir,$rside);
   my $both=0;
   my $rev=0;
   my @reas=();
   if ($q=~/(<=>|<=|=>)/){
      print "2-sided";
      ($lside,$dir,$rside)=split(/\s*(<=>|<=|=>)\s*/,$q);
      if ($dir eq "<=>"){
         $rev=1;
      } elsif ($dir eq "<="){
         _swap(*lside,*rside);
      }
   } else {
      print "1-sided";
      $lside=$q;
      $both=1;
   }
   $lside=[split(/\s*\|\s*/,$lside)];
   $rside=[split(/\s*\|\s*/,$rside)];
   for my $a ($lside,$rside) {for my $o (@$a){$o=[split(/\s+/,$o)]}}
   print dmpml({lside=>$lside,rside=>$rside});
   for my $r (@{$self->{reaID}}){
      my $s=$self->substratesID($r);
      my $p=$self->productsID($r);
      if ($both){
         push(@reas,$r) if _testcondset([@$s,@$p],$lside);
      } else {
         if ($rev){
            push(@reas,$r) if (_testcondset($s,$lside) && _testcondset($p,$rside)) || (_testcondset($p,$lside) && _testcondset($s,$rside));
         } else {
            push(@reas,$r) if _testcondset($s,$lside) && _testcondset($p,$rside);
         }
      }
   }
   return \@reas;
}
sub _testcondset{
   my $set=shift;
   my $conds=shift;
   my %seth=();
   @seth{@$set}=(1) x scalar(@$set);
   my $met=1;
   for my $c (@$conds){
      for my $subc (@$c){
         if ($subc=~/^\!/){
            my $rsubc=substr($subc,1);
            $met=0 if exists $seth{$rsubc};
         } else {
            $met=0 if ! exists $seth{$subc};
         }
      }
   }
#  print "test ".dmp($conds)." in",keys %seth,"ok\n" if $met;  
   return $met;
}

sub printReas{
   my $self=shift;
   my @reas=@_; # query string
   @reas=@{$reas[0]} if (ref($reas[0]) eq 'ARRAY');
#  @reas=@{$self->{reaID}} if ! defined $reas[0] ;
   my @arr=("<=","<=>","=>");
   for my $r (@reas){
      my $s=$self->substratesID($r);
      my $p=$self->productsID($r);
      my $ss=$self->subs_stoichID($r);
      my $ps=$self->prod_stoichID($r);
      for (my $i=0;$i<scalar(@$s);$i++){
         print " + " if $i>0;
         print $ss->[$i].$s->[$i];
      }
      print " ".$arr[$self->revID($r)+1]." ";
      for (my $i=0;$i<scalar(@$p);$i++){
         print " + " if $i>0;
         print $ps->[$i].$p->[$i];
      }
      print "\n";
   }
}

=head2 convert... (methods)

=over

=item $Net->convert[Compound|Reactions][ID|Idx|Msk]To[ID|Idx|Msk](array)

=item $Net->convertCompoundsIdxToMsk(array)

=item ...

=back

Convert Compounds or Reaction between various formats. ID represents an array-ref of compounds or reaction IDs. Idx is a list of compound or reaction indices in the network and Msk signifies a CArray containing a non-zero value at each index whose reaction/compound is present.

=head3 Example:

   $s=scope($Net->convertCompoundsIDtoMsk([C00002,C00001]));
   print @{$Net->convertCompoundsMsktoID($s);

=head3 short forms

=over

=item [c|r][ID|IX|MS][ID|IX|MS](array)

=back

where c or r indicates compounds or reactions and ID=ID, IX=Idx and MS=Msk.

   $s=scope($Net->cIDMS([C00002,C00001]));
   print @{$Net->cMSID($s);

=cut

#---compound converters
sub comIDs{
   return $_[0]->{comID};
}
*getCompoundIDs=*comIDs;
sub comID{
   return $_[0]->{comID}->[$_[1]];
}
*getCompoundID=*comID;
sub comIX{
   return $_[0]->{comhash}->{$_[1]};
}
*getCompoundIdx=*comIX;
sub cIDMS{
   my $self=shift;
   my $msk=new Scopes::ArrayU8($self->numC());
   $msk->selectPut(new Scopes::ArrayS([@{$self->{comhash}}{@{$_[0]}}]),1);
   return $msk;
}
*convertCompoundsIDToMsk=*cIDMS;
sub cIXMS{
   my $self=shift;
   my $msk=new Scopes::ArrayU8($self->numC());
   $msk->selectPut(new Scopes::ArrayS(shift()),1);
   return $msk;
}
*convertCompoundsIdxToMsk=*cIXMS;
sub cMSIX{
   my $self=shift;
   my $msk=shift;
   return $msk->which()->toa();
}
*convertCompoundsMskToIdx=*cMSIX;
sub cIDIX{
   my $self=shift;
   return [@{$self->{comhash}}{@{shift()}}];
}
*convertCompoundsIDToIdx=*cIDIX;
sub cMSID{
   my $self=shift;
   my $msk=shift;
   return [@{$self->{comID}}[@{$msk->which()->toa()}]];
}
*convertCompoundsMskToID=*cMSID;
sub cIXID{
   my $self=shift;
   return [@{$self->{comID}}[@{shift()}]];
}
*convertCompoundsIdxToID=*cIXID;
sub convertCompoundsMskToIDValHash {
   my $self=shift;
   my $msk=shift;
   my $options=shift;
   my $idx=0;
   my $ret={};
   foreach my $val (@{$msk->toa()}) {
       $idx++;
       next unless ($val || $options->{includezero});
       $ret->{$self->{comID}->[$idx-1]}=$val;
   }
   return $ret;
}
sub chIXID{ # channel indexes to ids
   my $self=shift;
   return [@{$self->{chnID}}[@{shift()}]];
}

#---reaction converters

sub reaIDs{
   return $_[0]->{reaID};
}
*getReactionIDs=*reaIDs;
sub reaID{
   return $_[0]->{reaID}->[$_[1]];
}
*getReactionID=*reaID;
sub reaIX{
   return $_[0]->{reahash}->{$_[1]};
}
*getReactionIdx=*comIX;
sub rIDMS{
   my $self=shift;
   my $msk=new Scopes::ArrayU8($self->numR());
   $msk->selectPut(new Scopes::ArrayS([@{$self->{reahash}}{@{$_[0]}}]),1);
   return $msk;
}
*convertReactionsIDToMsk=*rIDMS;
sub rIXMS{
   my $self=shift;
   my $msk=new Scopes::ArrayU8($self->numR());
   $msk->selectPut(new Scopes::ArrayS(shift()),1);
   return $msk;
}
*convertReactionsIdxToMsk=*rIXMS;
sub rMSIX{
   my $self=shift;
   my $msk=shift;
   return $msk->which()->toa();
}
*convertReactionsMskToIdx=*rMSIX;
sub rIDIX{
   my $self=shift;
   return [@{$self->{reahash}}{@{shift()}}];
}
*convertReactionsIDToIdx=*rIDIX;
sub rMSID{
   my $self=shift;
   my $msk=shift;
   return [@{$self->{reaID}}[@{$msk->which()->toa()}]];
}
*convertReactionsMskToID=*rMSID;
sub rIXID{
   my $self=shift;
   return [@{$self->{reaID}}[@{shift()}]];
}
*convertReactionsIdxToID=*rIXID;
sub convertReactionsMskToIDValHash {
   my $self=shift;
   my $msk=shift;
   my $options=shift;
   my $idx=0;
   my $ret={};
   foreach my $val (@{$msk->toa()}) {
       $idx++;
       next unless ($val || $options->{includezero});
       $ret->{$self->{reaID}->[$idx-1]}=$val;
   }
   return $ret;
}

#syntax: $net->addCompartments(name,[compound ID list])
sub addCompartment{
   my $self=shift;
   my $name=shift;
   $self->{compartments}->{$name}={};
   $self->addToCompartment($name,@_);
}
sub addToCompartment{
   my $self=shift;
   my $name=shift;
   my $coms=[@_];
   if (ref($coms->[0]) eq 'ARRAY'){ 
      $coms=$coms->[0];
   }
   for my $c (@$coms){
      if (exists $self->{comIDtoComp}->{$c}){
         delete $self->{compartments}->{$self->{comIDtoComp}->{$c}}->{$c};
      }
      $self->{compartments}->{$name}->{$c}=1;
      $self->{comIDtoComp}->{$c}=$name;
   }
}
sub getCompartment{
   my $self=shift;
   my $com=shift;
   return $self->{comIDtoComp}->{$com};
}
sub getCompartments{
   my $self=shift;
   my $name=shift;
   return [keys %{$self->{compartments}}];
}
sub getCompartmentComs{
   my $self=shift;
   my $name=shift;
   return [keys %{$self->{compartments}->{$name}}];
}
sub getCompartmentReas{
   my $self=shift;
   my $name=shift;
   my $mode=shift; # only for internal use of getCompartmentBoundReaList
   my @list=();
   for my $r (@{$self->reaIDs()}){
      my $allin=1;
      my $allout=1;
      my $ms=$self->metabolitesID($r);
      for my $c (@$ms){
         if ($self->{comIDtoComp}->{$c} ne $name){
            $allin=0;
         } else {
            $allout=0;
         }
      }
      if ($mode){
         push(@list,$r) if !($allin || $allout);
      } else {
         push(@list,$r) if $allin;
      }
   }
   return \@list;
}
sub getCompartmentBoundReas{
   my $self=shift;
   my $name=shift;
   return $self->getCompartmentReaList($name,1);
}
sub isComInCompartment{
   my $self=shift;
   my $c=shift;
   my $name=shift;
   return $self->{comIDtoComp}->{$c} eq $name;
}
sub isReaInCompartment{
   my $self=shift;
   my $r=shift;
   my $name=shift;
   my $allin=1;
   my $ms=$self->metabolitesID($r);
   for my $c (@$ms){
      $allin=0 if ($self->{comIDtoComp}->{$c} ne $name);
   }
   return $allin; 
}

sub setReaLabel{
   my $self=shift;
   my $id=shift;
   my $name=shift;
   $self->{rlabels}->{$id}=$name;
}
sub setComLabel{
   my $self=shift;
   my $id=shift;
   my $name=shift;
   $self->{clabels}->{$id}=$name;
}
sub getReaLabel{
   my $self=shift;
   my $id=shift;
   return $self->{rlabels}->{$id};
}
sub getComLabel{
   my $self=shift;
   my $id=shift;
   return $self->{clabels}->{$id};
}

#---------helpers-------------
sub _cIXAR{
   my $idxs=shift;
   die "_cIXAR expects array ref" unless ref($idxs);
   my $ar=new Scopes::ArrayS($idxs);
   $ar->DISOWN() if shift;
   return $ar;
}
sub _cIXARF{
   my $idxs=shift;
   die "_cIXAR expects array ref" unless ref($idxs);
   my $ar=new Scopes::ArrayF($idxs);
   $ar->DISOWN() if shift;
   return $ar;
}
sub _getcomar{
   my $self=shift;
   return _cIXAR($self->_getcomidx(shift()));
}
sub _getcomidx{
   my $self=shift;
   my $arid=shift;
   my $ref=1;
   return unless $arid;
   if (ref($arid) ne 'ARRAY'){
      $arid=[$arid];
      $ref=0;
   }
   my $aridx=[];
   for my $id (@$arid){
      my $idx;
      if (exists $self->{comhash}->{$id}){
         $idx=$self->{comhash}->{$id};
      } else {
         $idx=$self->{comhash}->{$id}=scalar(@{$self->{comID}});
         $self->{c_net}->resize($self->numR,$idx);
         push(@{$self->{comID}},$id);
         $self->{c_net}->setName($idx,$id);
      }
      push(@$aridx,$idx);
   }
   return $aridx if $ref;
   return $aridx->[0];
}
sub _getreaar{
   my $self=shift;
   return _cIXAR($self->_getreaidx(shift()));
}
sub _getreaidx{
   my $self=shift;
   my $arid=shift;
   return unless $arid;
   my $ref=1;
   if (ref($arid) ne 'ARRAY'){
      $arid=[$arid];
      $ref=0;
   }
   my $aridx=[];
   for my $id (@$arid){
      my $idx;
      if (exists $self->{reahash}->{$id}){
         $idx=$self->{reahash}->{$id};
      } else {
         $idx=$self->{reahash}->{$id}=scalar(@{$self->{reaID}});
         $self->{c_net}->resize($idx,$self->numC);
         push(@{$self->{reaID}},$id);
      }
      push(@$aridx,$idx);
   }
   return $aridx if $ref;
   return $aridx->[0];
}
sub _getchnar{
   my $self=shift;
   return _cIXAR($self->_getchnidx(shift()));
}
sub _getchnidx{
   my $self=shift;
   my $arid=shift;
   my $ref=1;
   return unless $arid;
   if (ref($arid) ne 'ARRAY'){
      $arid=[$arid];
      $ref=0;
   }
   my $aridx=[];
   for my $id (@$arid){
      my $idx;
      if (exists $self->{chnhash}->{$id}){
         $idx=$self->{chnhash}->{$id};
      } else {
         $idx=$self->{chnhash}->{$id}=scalar(@{$self->{chnID}});
         push(@{$self->{chnID}},$id);
      }
      push(@$aridx,$idx);
   }
   return $aridx if $ref;
   return $aridx->[0];
}
sub _swap(**){local *a=shift;local *b=shift;my $h=$a;$a=$b;$b=$h}

sub toString {
    my $self = shift;
    my ( $m, $n ) = $self->size;
    my @sizes   = ();
    my $sizesum = 0;
    my $text    = '';
    if ( ( $m < 55 ) && ( $n < 80 ) ) {
        for ( my $i = 0 ; $i < $m ; $i++ ) {
            my $max = 0;
            for ( my $j = 0 ; $j < $n ; $j++ ) {
                my $size = length( abs( $self->at( $i, $j ) ) );
                if ( $size > $max ) {
                    $max = $size;
                }
            }
            push( @sizes, $max );
            $sizesum += $max + 1;
        }
        $sizesum--;

        #      print @sizes,"\n";
        #      print " ";
        #      for (my $i=0;$i<$sizesum;$i++){
        #         print "-";
        #      }

        $text .= "\n";
        for ( my $j = 0 ; $j < $n ; $j++ ) {
            $text .= "|";
            for ( my $i = 0 ; $i < $m ; $i++ ) {
                if ( $i != 0 ) { $text .= " " }
                my $el = $self->at( $i, $j );
                my $elt=abs($el);
                if    ( $el < 0 ) { $text .= RED }
                elsif ( $el > 0 ) { $text .= BLUE }
                elsif ( InSet($self->catalystsIX($i),$j) ) { $elt="c";$text .= GREEN }
                elsif ( InSet($self->inhibitorsIX($i),$j) ) { $elt="i";$text .= RED }
                else              { $text .= WHITE }
                if ( InSet($self->inhibitedIX($i),$j) ) { $text .= REVERSE }
                if ( InSet($self->propagatorsIX($i),$j) ) { $text .= REVERSE}
                $text .= " " x ( $sizes[$i] - length( $elt ) );
                $text .= $elt;
                $text .= RESET;
            }
            $text .= "|  ";
            if ( $self->{comID}->[$j] ) {
                my $name = $self->{clabels}->{ $self->{comID}->[$j] };
                $text .= ($name ? sprintf( '%.*s', 120 - $sizesum - 10, $name ) : $self->{comID}->[$j] . " ");
            }
            $text .= "\n";
        }

        #      print " ";
        #      for (my $i=0;$i<$sizesum;$i++){
        #         print "-";
        #      }
        $text .= "\n";
        my $maxsize = 0;
        for ( my $j = 0 ; $j < $m ; $j++ ) {
            my $len = length( $self->{reaID}->[$j] );
            $maxsize = $len if ( $len > $maxsize );
        }
        for ( my $i = 0 ; $i < $maxsize ; $i++ ) {
            for ( my $j = 0 ; $j < $m ; $j++ ) {

                #          if ($j==0){
                #             print spacer($sizes[$j]-1);
                #          }else{
                #             print spacer($sizes[$j]);
                #          }
                $text .= " " x ( $sizes[$j] );
                $text .= RED if !$self->selectedIX($j);
                my $letter = ( split( //, $self->{reaID}->[$j] ) )[$i];
                $text .= defined $letter ? $letter : " ";
                $text .= RESET  if !$self->selectedIX($j);
            }
            if ( $i == 0 ) {
                $text .= "   " . RED . "1" . RESET . " substrates";
            }
            if ( $i == 1 ) {
                $text .= "   " . BLUE . "1" . RESET . " products";
            }
            if ( $i == 2 ) {
                $text .= "   " . GREEN . "c" . RESET . " catalysts";
            }
            if ( $i == 3 ) {
                $text .= "   " . RED . "i" . RESET . " inhibitors";
            }
            if ( $i == 4 ) {
                $text .= "   " . RED . REVERSE . "1" . RESET . " inh. subs.";
            }
            if ( $i == 5 ) {
                $text .= "   " . BLUE . REVERSE . "1" . RESET . " propagators";
            }
            $text .= "\n";
        }
    }
    else {
        $text .=
            "Scopes::Net: "
          . scalar( keys %{ $self->{reahash} } )
          . " reactions, "
          . scalar( keys %{ $self->{comhash} } )
          . " compounds";
    }
    return $text;
}
1;
